/***
 * 编程实现 A 为 100 × 100 随机矩阵时,用 QR 和 Cholesky 分解求 x 的程序。eigen最大支持50x50的矩阵
 */
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <ctime> //计时

#define MAT_SIZE 100

int main()
{
    /* 求解方程组 Ax=b */
    Eigen::MatrixXd randm = Eigen::MatrixXd::Random(MAT_SIZE, MAT_SIZE);
    Eigen::Matrix<double, MAT_SIZE, MAT_SIZE> mat_a = randm * randm.transpose();
    Eigen::Matrix<double, MAT_SIZE, 1> v_b = Eigen::MatrixXd::Random(MAT_SIZE, 1);
    std::cout << "mat_a:" << std::endl
              << mat_a << std::endl;
    std::cout << "v_b:" << std::endl
              << v_b << std::endl;

    // 使用最基础的逆矩阵求解
    std::clock_t start = std::clock();
    auto x = mat_a.inverse() * v_b;
    std::cout << "逆矩阵求解 x= " << x.transpose() << std::endl;
    std::clock_t end = std::clock();
    std::cout << "时间开销: " << (end - start) / (double)CLOCKS_PER_SEC * 1000 << "ms" << std::endl;

    // 使用QR分解求解方程组
    start = std::clock();
    auto xx = mat_a.colPivHouseholderQr().solve(v_b);
    // auto xx = mat_a.householderQr().solve(v_b);
    std::cout << "QR分解求解 x= " << xx.transpose() << std::endl;
    end = std::clock();
    std::cout << "时间开销: " << (end - start) / (double)CLOCKS_PER_SEC * 1000 << "ms" << std::endl;

    // 使用楚列斯基分解分解
    start = std::clock();
    auto xxx = mat_a.ldlt().solve(v_b);
    std::cout << "QR分解求解 x= " << xxx.transpose() << std::endl;
    end = std::clock();
    std::cout << "时间开销: " << (end - start) / (double)CLOCKS_PER_SEC * 1000 << "ms" << std::endl;

    return 0;
}